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Abstract — A sum-of-squares is a polynomial that can be ex- 
pressed as a sum of squares of other polynomials. Determining if 
a sum-of-squares decomposition exists for a given polynomial is 
equivalent to a linear matrix inequality feasibility problem. The 
computation required to solve the feasibility problem depends 
on the number of monomials used in the decomposition. The 
Newton polytope is a method to prune unnecessary monomials 
from the decomposition. This method requires the construction 
of a convex hull and this can be time consuming for polynomials 
with many terms. This paper presents a new algorithm for 
removing monomials based on a simple property of positive 
semidefinite matrices. It returns a set of monomials that is 
never larger than the set returned by the Newton polytope 
method and, for some polynomials, is a strictly smaller set. 
Moreover, the algorithm takes significantly less computation 
than the convex hull construction. This algorithm is then 
extended to a more general simplification method for sum-of- 
squares programming. 

I. Introduction 

A polynomial is a sum-of-squares (SOS) if it can be 
expressed as a sum of squares of other polynomials. There 
are close connections between SOS polynomials and positive 
semidefinite matrices [3], [2], [4], [13], [11], [7], [12]. For 
a given polynomial the search for an SOS decomposition is 
equivalent to a linear matrix inequality feasibility problem. 
It is also possible to formulate optimization problems with 
polynomial sum-of-squares constraints [11], [12]. There is 
freely available software that can be used to solve these SOS 
feasibility and optimization problems [14], [8], [1], [6]. Many 
nonlinear analysis problems, e.g. Lyapunov stability analysis, 
can be formulated within this optimization framework [11], 
[12], [19], [20]. 

Computational growth is a significant issue for these 
optimization problems. For example, consider the search for 
an SOS decomposition: given a polynomial p and a vector 
of monomials z, does there exist a matrix Q >; such 
that p = z T Qzl The computation required to solve the 
corresponding linear matrix inequality feasibility problem 
grows with the number of monomials in the vector z. The 
Newton polytope [15], [18] is a method to prune unnecessary 
monomials from the vector z. This method is implemented in 
SOSTOOLs [14]. One drawback is that this method requires 
the construction of a convex hull and this construction itself 
can be time consuming for polynomials with many terms. 

This paper presents an alternative monomial reduction 
method called the zero diagonal algorithm. This algorithm is 
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based on a simple property of positive semidefinite matrices: 
if the (i, i) diagonal entry of a positive semidefinite matrix 
is zero then the entire i th row and column must be zero. 
The zero diagonal algorithm simply searches for diagonal 
entries of Q that are constrained to be zero and then prunes 
the corresponding monomials. This algorithm can be imple- 
mented with very little computational cost using the Matlab 
find command. It is shown that final list of monomials 
returned by the zero diagonal algorithm is never larger than 
the pruned list obtained from the Newton polytope method. 
For some problems the zero diagonal algorithm returns a 
strictly smaller set of monomials. Results contained in this 
paper are similar to and preceded by those found in the prior 
work [9], [21]. 

The basic idea in the zero diagonal algorithm is then 
extended to a more general simplification method for sum- 
of-squares programs. The more general method also removes 
free variables that are implicitly constrained to be equal 
to zero. This can improve the numerical conditioning and 
reduce the computation time required to solve the SOS pro- 
gram. Both the zero diagonal elimination algorithm and the 
simplification procedure for SOS programs are implemented 
in SOSOPT [1]. 

II. SOS Polynomials 

N denotes the set of nonnegative integers, {0, 1, . . .}, and 
N n is the set of n-dimensional vectors with entries in N. 
For aeN", a monomial in variables {xi, . . . ,x n } is given 
by x a :— x^x^ 2 ■ ■ -x"". a is the degree vector associated 
with the monomial x a . The degree of a monomial is defined 
as degx" :— Yn=i a *- ^ polynomial is a finite linear 
combination of monomials: 



(1) 



where c a e 1, c a / 0, and A is a finite collection of vectors 
in N n . R[x] denotes the set of all polynomials in variables 
{xi, . . . ,x n } with real coefficients. Using the definition of 
deg for a monomial, the degree of p is defined as degp := 
max oe ^ [dega; Q ]. 

A polynomial p is a sum-of-squares (SOS) if there exist 
polynomials {fi}f =1 such that p = £™i fi- The set of SOS 
polynomials is a subset of M.[x] and is denoted by S[x]. If 
p is a sum-of-squares then p(x) > Vx £ M. n , However, 
non-negative polynomials are not necessarily SOS [16]. 

Define z as the column vector of all monomials in vari- 



ables {xi, . . . , x n } of degree < d: Q 

Z ■— [l, El, X2, Xn, X\ , HIj, . . . , 2&, . . . , X n ] T (2) 

There are ft -1 ) monomials in n variables of degree k. 
Thus z is a column vector of length l z := 2~2k=o ( fe+ fe 1 ) = 
If / is a polynomial in n variables with degree 
< d then by definition / is a finite linear combination of 
monomials of degree < d. Consequently, there exists a £ 
such that / = a T z. 

Two useful facts from [15] are: 

1) If p is a sum-of-squares then p must have even degree. 

2) If p is degree 2d (d 6 N) and p = YT=i fi then 
deg fi < d Vi 

The following theorem, introduced as the "Gram Matrix" 
method by [4], [13], connects SOS polynomials and positive 
semidefinite matrices. 

Theorem 1: Let p € M[x] be a polynomial of degree 2d 
and z be the Z z x 1 vector of monomials defined in Equation|2] 
Then p is a SOS if and only if there exists a symmetric matrix 
Q e R hxl * such that Q y and p = z T Qz. 

Proof: If p is a SOS, then there exists polynomials 

{fijili such that p = YT=i fi- B y fact 2 above ' de § h< d 
for all i. Thus, for each fi there exists a vector, dj € R z , 
such that /j = af z. Define the matrix A S R z * xm whose i" 1 
column is <ij and define Q := AA T y 0. Then p = z T Qz. 

(<=) Assume there exists Q = Q T e R lzXlz such that 
Q y and p = z T Qz. Define m := rank(Q). There exists 
a matrix Ael'» xm such that Q = AA T . Let a, denote the 
i th column of A and define the polynomials fi :— z T <n. By 
definition of f it p = z T (AA T )z = YZi fi- ■ 

Determining if an SOS decomposition exists for a given 
polynomial p is equivalent to a feasibility problem: 

Find QhO such that p = z T Qz (3) 

Q is constrained to be positive semi-definite and equating co- 
efficients of p and z T Qz imposes linear equality constraints 
on the entries of Q. Thus this is a linear matrix inequality 
(LMI) feasibility problem. There is software available to 
solve for SOS decompositions [14], [8], [1]. These toolboxes 
convert the SOS feasibility problem to an LMI problem. 
The LMI problem is then solved with a freely available 
LMI solver, e.g. Sedumi [17], and an SOS decomposition 
is constructed if a feasible solution is found. These software 
packages also solve SOS synthesis problems where some 
of the coefficients of the polynomial are treated as free 
variables to be computed as part of the optimization. These 
more general SOS optimization problems are discussed fur- 
ther in Section [V] Many analysis problems for polynomial 
dynamical systems can be posed within this SOS synthesis 
framework [11], [12], [19], [20]. 

'Any ordering of the monomials can be used to form z. In Equation f2] 
x a precedes x@ in the definition of z if Aegx a < degx' 3 OR dega; a = 
deg x@ and the first nonzero entry of a — fi is > 0. 



III. Newton Polytope 

As discussed in the previous section, the search for an SOS 
decomposition is equivalent to an LMI feasibility problem. 
One issue is that the computational complexity of this LMI 
feasibility problem grows with the dimension of the Gram 
matrix. For a polynomial of degree 2d in n variables there 
are, in general, l z — ( n ^ d ) monomials in z and the Gram 
matrix Q is l z x l z . l z grows rapidly with both the number 
of variables and the degree of the polynomial. However, any 
particular polynomial p may have an SOS decomposition 
with fewer monomials. The Newton Polytope [15], [18] is an 
algorithm to reduce the dimension l z by pruning unnecessary 
monomials from z. 

First, some terminology is provided regarding polytopes 
[10], [5]. For any set A C W 1 , convhull(A) denotes the 
convex hull of A. Let C C R" be a convex set. A point 
a G C is an extreme point if it does not belong to the 
relative interior of any segment [01,03] C C. In other 
words, if 3oi,02 6 C and < A < 1 such that a = 
A01 + (1 — A)«2 then oi = = o. A convex polytope (or 
simply polytope) is the convex hull of a non-empty, finite 
set {ai, . . . ,o p } C W l . The extreme points of a polytope 
are called the vertices. Let C be a polytope and let V be 
the (finite) set of vertices of C. Then C = convhull(V) and 
V is a minimal vertex representation of C. The polytope C 
may be equivalently described as an intersection of a finite 
collection of halfspaces, i.e. there exists a matrix H € R Arx ™ 
and a vector g eR N such that C = {a£l" : Ha < g}. 
This is a facet or half-space representation of C. 

The Newton Polytope (or cage) of a polynomial p = 
J2aeA c ^ xa is defined as C{p) := convhull(^l) C R" [15]. 
The reduced Newton polytope is \C(p) := {^a : a £ 
C(p)}. The following theorem from [15] is a key result for 
monomial reduction. |3 

Theorem 2: If p = J2iLi fi then the vertices of C(p) are 
vectors whose entries are even numbers and C(/,) C ^C(p). 

This theorem implies that any monomial x a appearing in 
the vector z of an SOS decomposition z T Qz must satisfy 
a E \C(jp) fi N". This forms the basis for the Newton 
polytope method for pruning monomials: Let p be a given 
polynomial of degree 2d in n variables with monomial 
degree vectors specified by the finite set A. First, create the 
l z x 1 vector z consisting of all monomials of degree < d in n 
variables. There are l z = ( n ~^ d ) monomials in this complete 
list. Second, compute a half-space representation {a £ R n : 
Ha < g} for the reduced Newton polytope ^C(p). Third, 
prune out any monomials in z that are not elements of 
\C{p). This algorithm is implemented in SOSTOOLs [14]. 
The third step amounts to checking each monomial in z 
to see if the corresponding degree vector satisfies the half- 
plane constraints Ha < g. This step is computationally 

2 A polynomial p is a form if all monomials have the same degree. The 
results in [15] are stated and proved for forms. A given polynomial can 
be converted to a form by adding a single dummy variable of appropriate 
degree to each monomial. The results in [15] apply to polynomials by this 
homogenization procedure. 



very fast. The second step requires computing a half-plane 
representation for the convex hull of hA. This can be done 
in Matlab, e.g. with convhulln. However, this step can 
be time-consuming when the polynomial has many terms 
(.4 has many elements). The next section provides an al- 
ternative implementation of the Newton Polytope algorithm 
that avoids constructing the half-space representation of the 
reduced Newton polytope. 



Example: Consider the following polynomial 

p — 3xf — 2x\x2 + 7xf — ix±X2 



1 



(4) 



p is a degree four polynomial in two variables. The list of 
all monomials in two variables with degree < 2 is: 



Z = [l X\ X2 x\ X\X2 X2] 



(5) 



The length of z is l z = 6. An SOS decomposition of a degree 
four polynomial would, in general, include all six of these 
monomials. The Newton Polytope can be used to prune some 
unnecessary monomials in this list. 

The set of monomial degree vectors for p is A := 
{[4], [2], [g], [i], [0] ) [0]}. These vectors are shown 
as circles in Figure Q] The Newton Polytope C(p) is the 
large triangle with vertices {[q] , [®] , [2]}- Figure [2] shows 
the degree vectors for the six monomials in z (circles) 
and the reduced Newton polytope (large triangle). The re- 
duced Newton polytope \C{f) is the triangle with vertices 



{[§], 



[ J]}. By Theorem|2] X1X2 and x\ can not appear 



in any SOS decomposition of p because [}] , [%] £ |C(/). 
These monomials can be pruned from z and the search for 
an SOS decomposition can be performed using only the four 
monomials in the reduced Newton polytope: 

iT 



Z = [1 Xi X2 x\\ ' 



(6) 



The length of the reduced vector z is l z = 4. The SOS 
feasibility problem with this reduced vector z (Equation [3]l 
is feasible. The following matrix is one feasible solution: 



Q = 



r 1 

7-2 

0-2 4 

0-1 



1 


-1 

3 



(7) 



p is SOS since p = z T Qz and Q >z 0. 

IV. Zero Diagonal Algorithm 

The zero diagonal algorithm searches for diagonal entries 
of the Gram matrix that are constrained to be zero and then 
prunes the associated monomials from z. The remainder of 
the section describes this algorithm in more detail. 

As mentioned in Section |IlJ equating the coefficients of p 
and z T Qz leads to linear equality constraints on the entries 
of Q. The structure of these equations plays an important 
role in the proposed algorithm. Let z be the l z x 1 vector 
of all monomials in n variables of degree < d (Equation |2). 
Define the corresponding set of degree vectors as M := 
{ai,...,ai r } C N™. z T Qz is a polynomial in x with 
coefficients that are linear functions of the entries of Q: 



z T Qz = 



m rn 

EE 

i=i 3=1 



ai+otj 



(8) 
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Fig. 1. Newton polytope (large triangle) and monomial degree vectors 
(circles) 
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Reduced Newton Polytope 
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Fig. 2. Reduced Newton polytope (large triangle) and degree vectors for 
all monomials of degree = 0, 1,2 (circles) 



The entries of z are not independent: it is possible that 



= ZkZi for some i,j,k,l € {1, . . . ,l z }. The unique 
degree vectors in Equation [8] are given by the set 

M + M := {a € N" : 3a it aj E M s.t. a = a t + aj} (9) 

The polynomial z T Qz can be rewritten as: 



z T Qz 



E E QwK 

a£M+M \(i,j)eS a 



where S a := : «i + aj — a}. Equating the 

coefficients of p and z T Qz yields the following linear 
equality constraints on the entries of Q: 

c a a e A 
at A 



E = { 



(ID 



There exists A £ M. lxl * and b G K ( such that these equality 
constraints are given by Aq = 60 where q := vec(Q) is the 
vector obtained by vertically stacking the columns of Q. The 
dimension I is equal to the number of elements of M + M. 

3 In addition to the equality constraints due to p = z T Qz there are also 
equality constraints due to the symmetry condition Q = Q T . Some solvers, 
e.g. Sedumi [17], internally handle these symmetry constraints. 



The zero diagonal algorithm is based on two lemmas. 



Lemma 1: If S^a; 
Lemma 2: If p = 



{(i, i)} then 







2a; G A 
2cti £ A 



(12) 



"Qz, Q y 0, and Qi 



then p 



z T Qz where z is the (l z — 1) x 1 vector obtained by deleting 
the i th element of z and Q >; is the (l z — 1) x (l z — 1) 
matrix obtained by deleting the i t/l row and column from Q. 

Lemma Q] follows from Equation [TT] 52 ft , = {(i,i)} 
means that x ai ■ x ai is the unique decomposition of x 2ai as a 
product of monomials in z. There is no other decomposition 
of x 2ai as a product of monomials in z. In this case, 
p = z T Qz places a direct constraint on Qn that must hold 
for all possible Gram matrices. 

Lemma [2] follows from a simple property of positive 
semidefinite matrices: If Q >z and Qi : i — then Qij — 
Q j{ = for j = l,...,l z . If Qi,i = then an SOS 
decomposition of p, if one exists, does not depend on the 
monomial Zi and Zi can be removed from z. 

The zero diagonal algorithm is given in Table U The sets 
Mfc denote the pruned list of monomial degree vectors at 
the k th iterate. The main step in the iteration is the search 
for equations that directly constrain a diagonal entry Qj j to 
be zero (Step 6). This step can be performed very fast since 
it can be implemented using the find command in Matlab. 
Based on Lemma|2l if Q^j = then the monomial Zi and the 
i th row and column of Q can be removed. This is equivalent 
to zeroing out the corresponding columns of A (Step 7). 
This implementation has the advantage that A and b do not 
need to be recomputed for each updated set M k . Zeroing 
out columns of A in Step 7 also means that new equations 
of the form Qi^ = may be uncovered during the next 
iteration. The iteration continues until no new zero diagonal 
entries of Q are discovered. The next theorem proves that if 
p is a SOS then the decomposition must be expressible using 
only monomials associated with the final set M& , ■ Moreover, 
Mf-f C 7}C(p) n N™, i.e. the list of monomials returned 
by the zero diagonal algorithm is never larger than the list 
obtained from the Newton polytope method. In fact, there are 
polynomials for which the zero diagonal algorithm returns a 
strictly smaller list of monomials than the Newton polytope. 
The second example below provides an instance of this fact. 

1. Given : A polynomial p = 5Da6.A c aX a . 

2. Initialization : Set k = and M := {a i } l i z =l C N" 

3. Form Aq = b: Construct the equality constraint data, A £ ]R !x ^ and 

b £ K', obtained by equating coefficients of p = z T Qz. 

4. Iteration : 

5. Set Z = 0, k := k + 1, and M k := M k -i 

6. Search Aq = b: If there is an equation of the form Qi i = 
then set Mj. := Mfc\{aj} and Z = Z U I where I are the 
entries of q corresponding to the i th row and column of Q. 

1. For each j E Z set the j th column of A equal to zero. 

8. Terminate if Z = otherwise return to step 5. 

9. Return : M k , A, b 

TABLE I 

Monomial Reduction using the Zero Diagonal Algorithm 



Theorem 3: The zero diagonal algorithm terminates in 
a finite number of steps, kf, and Mf., C hC(p) D N n . 
Moreover, if p = YT=i fi then <?(/*) n N ™ ^ M k f - 

Proof: Mo has l z elements. The algorithm terminates 
unless at least one point is removed from Mfc. Thus the 
algorithm must terminate after kf < l z + 1 steps. 

To show Mkj, C \C(p) Pi W 1 consider a vertex ctfj of 
convhull(M / t / ). If there exists u,v G convhul^Mfcj. ) such 
that 2oti = u + v then u = v = on. This follows from on — 
h(u + v) and the definition of a vertex. As a consequence, 
S2a z = {(hi)}- By Lemma[T] 



C2ai 2cti G A 
2a>i<£A 



(13) 



Qi i ^ since a* was not removed at step 6 during the final 
iteration and thus 2a; G A C C(p). This implies that Oj G 
hC(p), i.e. \C(p) contains all vertices of convhul^Mfej.). 
Hence M kf C convhull(A/ fc/ ) C |C(p). 

Finally it is shown that n N" C M kf . C(/j) C 

\C{p) by Theorem |2] and \C{p) C convhull(M ) by the 
choice of M . Thus C(fA n N n C M . Let z be the vector 
of monomials associated with Mo. If p = fi tnen there 

exists a Q >z such that p = z T Qz. If the iteration removes 
no degree vectors then M^ f — Mo and the proof is complete. 
Assume the iteration removes at least one degree vector and 
let oti be the first removed degree vector. Based on Step 6, 
p = z T Qz constrains Q^j = 0. By Lemma [2] the monomial 
Zi cannot appear in any Hence C(fi) n N n C Mo\{cti}. 
Induction can be used to show C(/j) flN™ C M k holds after 
each step k including the final step kf. ■ 

This algorithm is currently implemented in SOSOPT [1]. 
The results in Theorem [3] still hold if Mo C W 1 is chosen to 
be any set satisfying \C[p) D N n C M . Simple heuristics 
can be used to obtain an initial set of monomials Mo with 
fewer than l z elements. Mo can then be used to initialize the 
zero diagonal algorithm. The next step is to construct the 
matrix A and vector b obtained by equating the coefficients 
of p and z T Qz. This step is required to formulate the LMI 
feasibility problem and it is not an additional computational 
cost associated with the zero diagonal algorithm. M kf con- 
tains the final reduced set of monomial degree vectors. If at 
least one degree vector was pruned then the returned matrix 
A and vector b may contain entire columns or rows of zeros. 
These rows/and columns can be deleted prior to passing the 
data to a a semi-definite programming solver. The next two 
examples demonstrate the basic ideas of the algorithm. 

Example: Consider again the polynomial in Equation |4] 
The full list of all monomials in two variables with degree 
< 2 consists of six monomials (Equation |5). Equating 
the coefficients of p and z T Qz yields the following linear 



equality constraints on the entries of Q: 

Qi.i + Q 
Qe,4 + Q 
Qe.i + Q 

Q$,2 + Q2,5 + Q 



Q2,l 


+ 


Ql,2 


= o, 


Qi,2 


+ 


Q2,4 


= o, 


Qa.i 


+ 


Ql,3 


= 0, 


Q5,4 


+ 


Q4,5 


= 0, 


Qe,i 


+ 


<?3,6 


= 0, 






= 0, 






Qi,i 


= 1 






Qc>,6 


= 



Q(>,2 + Q 2 ,6 + Q 5 ,3 + 

Qs,i + Qi,5 + 1 



?2,2 = 7 

?5,5 = 

?3,3 = 4 

?3,4 = -2 

?3,5 = 

>3,2 + Q2,3 = -4 

Q4.4 = 3 



1,4 + 
4,6 + 
1,6 + 
4,3 + 



A matrix A and vector b can be constructed to represent 
these equations in the form Aq = b. Note that Q$,q = 
and this implies that Q^ g = Qq.i = i = 1, . . . , 6 for any 
SOS decomposition of p. Thus the monomial zq = x\ can 
not appear in any SOS decomposition and it can be removed 
from the list. After eliminating x\ and removing the 6 th row 
and column of Q, the equality constraints reduce to: 



?2,1 +Ql,2 = 0, 
?4,2 + Q2,4 = 0, 
?3,1 + Ql,3 = 0, 
?5,4 + Q4,5 = 0, 
?5,3 + Q 3 ,5 = 
Ql.l = 1 



?4,1 + ' 



.4 + 



?2,2 = 7 
?5,5 = 



Q 3 ,3 = 4 
?5,2 + Q 2 ,5 + Q4,3 + Q3,4 = -2 
?5,1 + Ql,5 + Q:i,2 + Q 2 ,3 = -4 

Q 4 4 = 3 



Removing the 6 th row and column of Q is equivalent to 
zeroing out the appropriate columns of the matrix A. This 
uncovers the new constraint Q5 5 = which implies that 
the monomial £5 = X1X2 can be pruned from the list. 
After eliminating X1X2, the procedure can be repeated once 
again after removing the 5 th row and column of Q. No new 
diagonal entries of Q are constrained to be zero and hence 
no additional monomials can be pruned from z. The final list 
of monomials consists of four monomials. 

z = [1 xi x 2 x\] T (14) 

The Newton polytope method returned the same list. 

Example: Consider the polynomial p = x\ + x\ + x\x\. 
The Newton polytope is C(p) = convhull({[2] , [°], [|]}) 



The reduced Newton polytope 
convhull({[i], [?], [§]}). The 
corresponding to \C(p) nN™ is: 



2J ) 

is \C{p) 
monomial vector 



z := [xi, x 2 , x x x 2 , xfxl] 



(15) 



There are l z — 15 monomials in two variables with degree 
< 4. For simplicity, assume the zero diagonal algorithm is 
initialized with Mo := \C(p) n N n . Equating coefficients 
of p and z T Qz yields the constraint Q 3 3 = in the first 
iteration of the zero diagonal algorithm. The monomial z 3 = 
X1X2 is pruned and no additional monomials are removed at 
the next iteration. The zero diagonal algorithm returns M 2 = 
{[q] > [?] > [§]}• M 2 is a proper subset of \C{p) n N n . 
The same set of monomials is returned by the zero diagonal 
algorithm after 13 steps if Mq is initialized with the l z = 15 
degree vectors corresponding to all possible monomials in 
two variables with degree < 4. This example demonstrates 
that the zero diagonal algorithm can return a strictly smaller 
set of monomials than the Newton polytope method. 



V. Simplification Method for SOS Programs 

This section describes a simplification method for SOS 
programs that is based on the zero diagonal algorithm. A 
sum-of-squares program is an optimization problem with 
a linear cost and affine SOS constraints on the decision 
variables [14]: 



min c T u 



(16) 



subject to: dk(x, u) G E[x], k = 1, . . . N 



u G W are decision variables. The polynomials {ak}^ =1 are 
given problem data and are affine in it: 

a k {x, u) := a kfi (x) + a k ,i(x)ui H h a k , r (x)u r (17) 

Theorem [T] is used to convert an SOS program into a 
semidefinite program (SDP). The constraint a k (x, u) G S[x] 
can be equivalently written as: 

ak,o( x ) + a k,i(x)ui H h a kir (x)u r = z^Q k z k (18) 

Qk h (19) 

If max u [degafc(x,?/)] = 2c? then, in general, z k must 
contain all monomials in n variables of degree < d. Q k 
is a new matrix of decision variables that is introduced 
when converting an SOS constraint to an LMI constraint. 
Equating the coefficients of z^Q k z k and a k (x,u) imposes 
linear equality constraints on the decision variables u and 
Q k . There exists a matrix A G R ixm and vector b G M 1 such 
that the linear equations for all SOS constraints are given by 
Ay = b where 

V ■= [u T , vec(Q 1 ) T , vec(Q N ) T ] T 



(20) 



vec(Qk) denotes the vector obtained by vertically stacking 
the columns of Q k . The dimension m is equal to r + 
J2 k =i m k wnere Qk is m k x m k (k — 1,...,N). After 
introducing a Gram matrix for each constraint the SOS 
program can be expressed as: 

min c T u (21) 

subject to: Ay = b 
QkhO: k = l,...N 

Equation [21] is an SDP expressed in Sedumi [17] primal 
form, it is a vector of free decision variables and {Qk} k =i 
contain decision variables that are constrained to lie in the 
positive semi-definite cone. Sedumi internally handles the 
symmetry constraints implied by Q k = Q\- 

The SOS simplification procedure is a generalization of 
the zero diagonal algorithm. It prunes the list of monomials 
used in each SOS constraint. It also attempts to remove free 
decision variables that are implicitly constrained to be zero. 
Specifically, the constraints in some SOS programs imply 
both Ui > and u\ < 0, i.e. there is an implicit constraint 
that Ui = for some i. Appendix A.l of [19] provides 
some simple examples of how these implicit constraints 
can arise in nonlinear analysis problems. For these simple 
examples it is possible to discover the implicit constraints by 



examination. For larger, more complicated analysis problems 
it can be difficult to detect that implicit constraints exist. The 
SOS simplification procedure described below automatically 
uncovers some classes of implicit constraints u, = and 
removes these decision variables from the optimization. 
This is important because implicit constraints can cause 
numerical issues for SDP solvers. A significant reduction in 
computation time and improvement in numerical accuracy 
has been observed when implicitly constrained variables are 
removed prior to calling Sedumi. 

The general SOS simplification procedure is shown in 
Table HIl To ease the notation the algorithm is only shown for 
the case of one SOS constraint (N — 1). The extension to 
SOS programs with multiple constraints (N > 1) is straight- 
forward. The algorithm is initialized with a finite set of vec- 
tors Mq C N". The Newton polytope of a(x, u) depends on 
the choice of u so Mo must be chosen so that it contains all 
possible reduced Newton polytopes. One choice is to initial- 
ize Mq corresponding to the degree vectors of all monomials 
in n variables and degree < 2d := max u [degafc(a;, u)]. A 
and b need to be computed when formulating the SDP so 
this step is not additional computation associated with the 
simplification procedure. The last pre-processing step is the 
initialization of the sign vector s. The entries of Sj are +1, 
— 1, or if it can be determined from the constraints that 
Hi is > 0, < or = 0, respectively. Sj =NaN if no sign 
information can be determined for yi. If y\ corresponds to a 
diagonal entry of Q then s, can be initialized to +1. 

The main iteration step is the search for equations that 
directly constrain any decision variable to be zero (Step 
7a). This is similar to the zero diagonal algorithm. The 
iteration also attempts to determine sign information about 
the decision variables. Steps 7b-7d update the sign vector 
based on equality constraints involving a single decision 
variable. For example, a decision variable must be zero if the 
decision variable has been previously determined to be < 
and the current equality constraint implies that it must be > 
(Step 7c). These decision variables can be removed from the 
optimization. Step 8 processes equality constraints involving 
two decision variables. The logic for this case is omitted due 
to space constraints. The processing of equality constraints 
can be performed very fast using the find command in 
Matlab. Steps 9 and 10 prune monomials and zero out 
appropriate columns of A. The iteration continues until no 
additional information can be determined about the sign of 
the decision variables. 

VI. Conclusions 

The Newton polytope is a method to prune unnecessary 
monomials from an SOS decomposition. The method re- 
quires the construction of a convex hull and this can be 
time time consuming for polynomials with many terms. 
This paper presented a zero diagonal algorithm for pruning 
monomials. The algorithm is based on a simple property of 
positive semidefinite matrices. The algorithm is fast since it 
only requires searching a set of linear equality constraints 
for those having certain properties. Moreover, the set of 



1. Given: Polynomials in variables x. Define 

a(x, u) := ao(x) + a\ (x)u\ + ■ ■ ■ + a r (x)u r 

2. Initialization: Set k = and choose a finite set Mo := {oiKLi 

C N n such that [u ueH r \C{a{x, u))] nN"CM . 

3. Form Ay = b: Construct the equality constraint data, A £ R ix ( r + m ) 

and b £ K 1 obtained by equating coefficients of a(x, u) = z T Qz 
where z := [x ai , . . . , x am ] and y := [u T , vec(Q) T ] T . 

4. Sign Data: Initialize the I X 1 vector s to be Sj = +1 if yi 

corresponds to a diagonal entry of Q. Otherwise set s, = NaN. 

5. Iteration : 

6. Set Z = 0, S = 0, k := k + 1, and M k := M k _ 1 

7. Process equality equality constraints of the form aijyj = bi 
where a^j ^ 

7a. If bi = then set sj = and Z = Z U j 

7b. Else if Sj =NaN then set Sj = sign(aijbj) and S = S U j 

7c. Else if Sj = —1 and sign(ciijbi) = +1 then set Sj = 

and S = S U j 

7d. Else if Sj = +1 and sign(a^ jbi) = — 1 then set Sj = 

and S = S U j 

8. Process equality equality constraints of the form 

a i>ji2/?i + a i,32yj2 = bi- 

9. If for any j £ Z, yj corresponds to a diagonal entry Qi^i 
then set M k := M k \{ai} and Z = Z U X where X are the 
entries of y corresponding to the i th row and column of Q. 

10. For each j £ Z set the j th column of A equal to zero. 

1 1 . Terminate if Z = and S = otherwise return to step 6. 

12. Return : M k , A, b, s 

TABLE II 

Simplification Method for SOS Programs with One 
Constraint 

monomials returned by the algorithm is a subset of the set 
returned by the Newton polytope method. The zero diagonal 
algorithm was extended to a more general reduction method 
for sum-of-squares programming. 
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